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Abstract 


Galaxies throughout the last 12 Gyr of cosmic time follow a sin- 
gle, universal relation that connects their star-formation rates 
(SFRs), stellar masses (M,) and chemical abundances [1, 2]. 
Deviation from these fundamental scaling relations would imply 
a drastic change in the processes that regulate galaxy evolution. 
Observations have hinted at the possibility that this relation 
may be broken in the very early universe [3, 4]. However, until 
recently, chemical abundances of galaxies could be only mea- 
sured reliably as far back as redshift z — 3.3 [5]. With JWST, 
we can now characterize the SFR, M,, and chemical abundance 
of galaxies during the first few hundred million years after the 
Big Bang, at redshifts z — 7 — 10. Here we show that galaxies 
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at this epoch follow unique SFR-M,—main-sequence and mass- 
metallicity scaling relations, but their chemical abundance is 
a factor of three lower than expected from the fundamental- 
metallicity relation of later galaxies. These findings suggest 
that galaxies at this time are still intimately connected with 
the intergalactic medium and subject to continuous infall of 
pristine gas which effectively dilutes their metal abundances. 


Main text 


Chemical enrichment is a key ingredient in the processes driving galaxy for- 
mation and evolution, but is also a ledger containing information on the early 
star-formation histories and stellar populations of the first galaxies [6, 7]. 
We identified galaxies at z > 7, ie. within the first 750 Myr after the Big 
Bang, through a blind and uniform selection of targets from recently-obtained 
public datasets with JWST Near-Infrared Spectrograph (NIRSpec) prism spec- 
troscopy. These include the gravitational lensing clusters Abell 2744 (JWST 
DD-2756, PI: Chen) and RXJ-2129 (JWST DD-2767, PI: Kelly), and the Cos- 
mic Evolution Early Release Science (CEERS) survey (JWST ERS-1345, PI: 
Finkelstein). We included all galaxies for which at least the nebular emission 
lines [O 111] \A4959, 5007 were detected, requiring S/N> 3 per wavelength ele- 
ment in the spectral regions covering the lines. Such a selection, in combination 
with the photometric identification of the targeted galaxies, has been shown to 
sample the intrinsic ultraviolet luminosity function of Lyman-break galaxies at 
z = 5-7, and is therefore likely representative of the underlying high-redshift 
galaxy population [8]. The spectroscopic and photometric data covering these 
fields were obtained from the MAST archive and reduced using grizli and 
custom-made pipelines (see Methods). We identify 16 targets matching the 
selection criteria, as summarized in Table 1. An example of the spectroscopic 
and imaging data obtained for each source is presented in Fig. 1. Throughout 
the paper we assume concordance ACDM cosmology [9]. 

We scale the optimally-extracted 1D spectra to match the photometry for 
each source using a wavelength-dependent polynomial function (see Methods). 
This is to improve the absolute flux calibration of the spectra and correct for 
potential slit-losses. We determine the spectroscopic redshift of each galaxy 
by fitting a continuum model with Gaussian line profiles of the most promi- 
nent nebular emission lines imposed on to each spectrum. The redshifts span 
Zspec = 7.11 — 9.50, corresponding to 500 — 750 Myr after the Big Bang, as 
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summarized in Table 1. We measure line fluxes for each transition modelled 
with the Gaussian line profiles in the photometrically-calibrated spectra. In 
all cases, we detect the [O11] 444959, 5007 doublet (by selection) and in the 
majority of cases we detect the Balmer line H8 and the [O r1] 443726, 3729 
doublet (14 and 12 out of the 16 galaxies, respectively). The [O 1] doublet is, 
however, not resolved in our spectra so we therefore only consider the sum 
of the doublet transitions, hence denoted [O 11] 43727. We are also not able 
to resolve the auroral [O 11] 44363 line transition from Hy, which is also only 
marginally detected in a few cases, except for the galaxy at z — 9.501 (see 
Methods). This feature will generally be resolved with the NIRSpec prism 
only for the highest redshift sources, where the line transition is located in the 
redder, higher wavelength-resolution region. 

We model the spectral energy distribution (SED) of each galaxy using the 
Bayesian Analysis of Galaxies for Physical Inference and Parameter EStima- 
tion (BAGPIPES) software package [10]. We infer the physical properties by 
jointly modeling the spectroscopic and photometric data, covering 0.7 — 5.2 um 
(~ 1000 — 5000A in the rest-frame) for each source. We assume a non- 
parametric star-formation history (SFH), since constant SFHs will typically 
underestimate the stellar mass build-up from older populations [11]. Indeed, 
spatially resolved galaxies at z > 7 have been observed to host older stellar 
populations that are otherwise outshone in integrated analyses by the younger, 
burstier stellar population [12]. We correct the SED outputs for the magnifi- 
cation factors due to lensing listed in Table 1, assuming the GLAFIC-v4 lens 
model for the Abell 2744 cluster, and a specific lens model determined for 
the RXJ-2129 galaxy cluster (see Methods). We infer stellar masses in the 
range M, = 107? — 1019? Mo, mass-weighted ages of the stellar populations of 
Tmass = 200—400 Myr, and low dust attenuation (typically Ay < 0.1 mag). The 
stellar masses and mass-weighted ages are consistent with previous estimates of 
galaxies at z > 7 assuming a similar flexible SFH, but are substantially larger 
(0.5-1dex) than inferred from models assuming a likely too simple constant 
SFH (see Methods). We apply an additional 0.2 dex systematic uncertainty to 
the statistical errors on the derived stellar masses in the following analysis, 
due to the systematic uncertainties stemming from the exact choice of SED 
modelling technique used [13]. 

We infer the star-formation rate (SFR) for each galaxy based on the HG 


flux measurements [14] as 


SFRug(Mo/yr) = 5.5 x 10" Lyg(erg/s) x aana ; (1) 
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assuming a Kroupa initial mass function (IMF) [15] and the theoretical 
fnajng = 2.86 ratio from the Case B recombination model at T, = 10^ K [16]. 
We measure SFRs in the range 2— 100 M; yr^!, based on the photometrically- 
calibrated spectra and taking into account the relevant magnification factors. 
The statistical uncertainties from the H£ line flux measurements are typically 
10— 3096, whereas the uncertainties from the choice of the IMF is = 20 — 30%. 
We therefore conservatively assume a 0.4 dex total uncertainty for each SFR 
measurement as summarized in Table 1. Combined with the stellar masses 
inferred from the photometry, we find that the identified galaxies at z — 7 — 10 
all seem to follow the same SFR-M, star-forming galaxy “main-sequence” 
(SFMS), quantified by log(SFR/Mo yr ^!) = 0.7 x log(M,/Mo) — 5.2 with a 
scatter of 0.4 dex (see Fig. 2). This is substantially higher, by 50, than the 
observed local, z ~ 0 SFMS [17], but consistent with previous mass-complete 
estimates of the SFMS at z ~ 6 [18]. 

We compare our observational results to predictions of the SFMS at 
z = 7 — 10 from recent cosmological simulations in Fig. 2, including: Astraeus 
[19], the First Light and Reionisation Epoch Simulations (FLARES) [20], First- 
Light [21], and the TNG50 simulations [22]. The simulations appear to be in 
good agreement internally, and are overall consistent with our data. We find 
marginal evidence for a flatter slope of the SFMS at z = 7 — 10 than predicted 
by simulations. However, this is mainly based on the two least massive galax- 
ies in our sample. The simulations are able to reproduce our results for the 
majority the sample, at stellar masses > 10° Mo. 

To determine the gas-phase metallicity through the inferred oxygen abun- 
dance, 12 + log(O/H), we rely on the strong-line diagnostics of the most 
prominent nebular emission lines as we are, in most cases, not able to resolve or 
securely constrain the auroral [O 111] 4363 line transition. We adopt the strong- 
line calibration based on the [O 11] A5008/H8 line ratio (denoted R3), con- 
strained over a large range of gas-phase metallicities, 12--log (O/H) — 6.9 — 8.9 
[23]. This particular line ratio was chosen due to the proximity of [O 111] A5008 
and Hf in wavelength space, such that the derived line ratios are not affected 
by systematic issues due to the flux calibration of the spectra or dust attenua- 
tion. Further, it has been established with a small empirical scatter (0.16 dex), 
and is supported by recent theoretical considerations [23] and empirical con- 
straints based on the detection of the auroral [O 11] M363 line in a sizable 
sample of galaxies at z > 5 [24]. Specifically, we adopt the “Large EW” R3 
calibration based on the observed equivalent width of H8, EWy, > 200 À, of 


the galaxies in our sample, representing the stronger ionization states. For the 
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galaxy RXJ-z9500, we place lower bounds on the inferred metallicity based 
on the marginal detection of [O 111] M363 using the direct T, method (Meth- 
ods), consistent with that inferred from the R3 diagnostic. For a subset (3 out 
of 16) of the galaxies we infer the metallicities using other calibrations, since 
their [O 111] A5008/H line ratios were at the turnover of the R3 calibration and 
therefore do not allow us to robustly constrain the metallicity (see Methods). 

'The mass-metallicity relation of the galaxies at z — 7 — 10 are shown in 
Fig. 3. We observe that they all follow the same approximate log-linear rela- 
tion, quantified as 12 + log(O/H) = 0.33 x log( M, / Mc) + 4.85, with a 0.3 dex 
scatter. The observed galaxies show substantially lower gas-phase metallici- 
ties, at 50 confidence, than local star-forming galaxies at z ~ 0 at equivalent 
stellar masses [2]. We again compare our observations to predictions of the 
mass-metallicity relation of galaxies at z — 7 — 10 from simulations, includ- 
ing constraints from Astraeus [25], FirstLight [26], FLARES [27], and TNG50 
[22]. Overall, the simulations show a substantial scatter in their predictions 
of the mass-metallicity relation at z — 7 — 10. Our observations are able to 
put strong constraints on the stellar and chemical build-up in these various 
simulations, and seem to best match the prescriptions used in Astraeus and 
FirstLight. The slopes of the empirical relation and those predicted from sim- 
ulations are generally consistent, except for FLARES which predicts a more 
rapid chemical enrichment with increasing stellar mass. 

With this unique dataset, we can for the first time test whether high- 
redshift galaxies follow the same fundamental mass-metallicity-SFR relation 
(FMR) observed for galaxies within the last 12 Gyr of cosmic time, i.e. from 
redshifts z = 0 — 3 [1, 2, 5]. Fig. 4 shows the metallicity offsets from the 
local FMR based on the inferred SFR and M, for each of the galaxies at 
z = 7 — 10 in our sample. We observe systematically lower metallicities (Zops) 
than predicted from the FMR (Zpyn), with a weighted mean and standard 
deviation of Zobs — ZFMR = —0.50+0.05 at z = 7— 10. The high-redshift FMR 


recovered from the Astraeus simulation [25] is generally consistent with our 


observations and predict an even larger offset from the local FMR for galaxies 
at z = 7 — 10 at low stellar masses or comparably higher SFRs (see Fig. 4). 
There have already been some early results hinting that the local FMR relation 
may not hold at redshifts z = 3.5 [1, 3, 4]. Most of these observations relied 
on limited sample sizes, however, and the majority of sources were starburst 
galaxies, making them overall unrepresentative of the typical galaxy population 
at these redshifts. The first JWST study of the three galaxies at z = 7 observed 
as part of the ERO was also inconclusive [28]. The results presented here are 
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thus the first to robustly verify this offset and the divergence from the local 
FMR at early cosmic times. 

'The drastic drop in metallicity in these early universe galaxies could be 
explained by significant accretion from the intergalactic medium that rapidly 
replenishes the galaxies with continuous infall of neutral, pristine gas. This 
process effectively dilutes the metals, reducing the overall chemical abundance 
[7]. Gas accretion from the intergalactic medium is expected theoretically to be 
an important process for feeding star-forming galaxies in the first ~ 1 — 2 Gyr 
of cosmic time [29]. Recent observations of abundant gas reservoirs in galaxies 
at z > 6 [30] support this hypothesis. In this scenario, it is the exhaustion of 
these pristine, intergalactic gas reservoirs that is responsible for the peak and 
rapid drop in the cosmic star-formation rate observed at redshift z ~ 1-3 [31]. 
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Table 1: Physical properties of the primary sample galaxies at z = 7 — 10. 


Galaxy ID Zspec SFRuyg log M, Ay Tmass 12 + log(O/H) m 


(Mo yr) (Mo) (mag) (Myr) 
RXJ-z9500 9.5008 33.91353 8.841902 0.027005 221+13  7.56t916 (R3) 99 
RXJ-z8149 8.1496 — 7.67548 — 819*005  0.05*009 262+17  7.29+0:22 (R3) 2.25 
RXJ-z8152 8.1523  124*155  8.72+0:02  0.04*005 266719 — T.68'019 (R3) 14640. 
Abell-z7878 7.8783 < 9.8 9.054906 0.2RTU detur > 7.30 (R3) 133-04 
Abell-27885 7.8854  8.8+13:3 Bee uP eS  o.10t002 265714 9770 25 (O32) 2120. 
CEERS-z7105 7.1051 5.1737 8.454 602 0.02+9:02 30412? ^ "TS UIS (R3) -— 
CEERSz7179 7.1702  7.8* Af. 833001 0.06* 0-05 20813) T4974 ; R3) => 
CEERS-27175 7.1751 14.344? 8.66102 0 021001 30579 742 91 R3) =s 
CEERS-z7167 7.1683  3.0!1?  840t003 0.031902  296*72  7.86+0:23 (R3) -- 
CEERS-z8684 8.6849 88.411336 10.0'001 0.51799; 398'2  — 8.06'017 (R3) -— 
CEERS-z7789 7.7805 «dd Shor eae DE ee OTI > 8.06 (R3) -- 
+28.7 +0:03 +0.01 +16 +0.21 
CEERS-z8718 8.7182 19.0777 — 9.05*002  0.02*0 0:  242*19 — 7.55 10-27 (R3) =r 
CEERS-z7832 7.8828 —17.3!*104 — 9.07700. o10t002 — 44779 — 7.75017 (R3) -— 
CEERS-z8612 8.6142 414 9 VET 0.00 o'p 234726 8 007034 (032) — 
CEERS-8172 81724 1.97779  904*9 WR, 261715 7.62 03s (R2) 25 
CEERS-z7453 7.4536  3.0'1$  — 8.64t00 0.021002 9841! 7.427055 (R3) -— 


Notes. The listed properties are not corrected for the magnification factor provided in Col. 8. 
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Fig. 1: Photometric and spectroscopic data of CEERS-z7382. Panel (a): False- 
color JWST/NIRCam RGB image zoomed in on the example galaxy (blue: 
F150W; green: F277W; red: F444W). The image scale and corresponding phys- 
ical size at z = 7.8328 is marked as well. Panel (b): Full NIRSpec/prism 
spectrum covering 0.7um to 5.2um (cyan) and associated error spectrum 
(grey). Panel (c): Zoom-in on the spectral region covering the nebular emis- 
sion lines from the [O 111] AA4960, 5008 doublet and H£. The local best-fit line 
and continuum model is shown by the black curve. 


Springer Nature 2021 ATEX template 


BN z=7-10 NIRSpec/prism 
(This work) 


FLARES 
TNG-50 
FirstLight 
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Star-formation rate, log(SFR/M o vr") 


7.0 FS 8.0 8.5 9.0 9.5 10.0 10.5 
Stellar mass, log(M , /M o) 


Fig. 2: The star-forming galaxy SFR-M, main-sequence at z — 7 — 10. Red 
squares show the primary sources analyzed in this work where the errorbars 
represent the combined statistical and systematic uncertainties. The best-fit 
linear relation is shown by the black curve, with the empirical scatter repre- 
sented by the grey-shaded region. The colored regions represent the predicted 
SFR-M, relations from z — 7 to z — 10, extracted from the simulations 
denoted in the legend. The bottom dark-grey curve shows the star-forming 
galaxy SFR-M, main-sequence at z — 0 [17]. 
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Fig. 3: The stellar mass-metallicity relation of galaxies at z — 7 — 10. The 
symbol notation for our observations and the simulations are identical to Fig. 2, 
except for the dark-grey dashed and continuous line which represents the mass- 
metallicity relation of galaxies at z = 0 from [2] (extrapolated below M, = 
108 Mo). 
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Fig. 4: Top panel: The fundamental-metallicity relation (FMR) of galax- 
ies at z = 0 — 3 [2, converted to a Kroupa IMF], extrapolated below 
log(M,/Mo) — 0.55 x log(SFR/Mo yr!) = 8.5. The symbol notation follows 
Fig. 2, with the functional form of the high-redshift FMR predicted by the 
Astraeus simulations [25] shown in blue. Bottom panel: Offset from the local 


FMR. The black solid line and grey-shaded region represent 
mean and standard deviation of the offset (A log(Z/Zo) = 0.504 
primary sample of galaxies at z = 7 — 10. 


the weighted 
- 0.05) for the 
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Methods 


Cosmology. Throughout this work, we assume concordance, flat ACDM cos- 
mology with Qm = 0.315, Qa = 0.685, and Hp = 67.8kms ! Mpc^! [9]. We 
use the cosmology distance calculator from astropy [32] to infer the luminosity 
distances Dr, to the given redshifts. 


Observations and data reduction. This work is primarily based on the 
photometric and spectroscopic data taken as part of the public UNCOVER 
(GO-2561), GLASS (ERS-1324), and CEERS (ERS-1345) surveys, and includ- 
ing observations obtained through the JWST Director's Discretionary time 
(DD-2756, PI: Chen and DD-2767, PI: Kelly). The NIRCam observations were 
typically taken in six broad-band filters (F115W, F150W, F200W, F277W, 
F356W, F444W). The photometric data are extracted from the catalog pre- 
sented by Brammer et al. (in prep), which provides a compilation of the JWST 
ERS photometric data released to date. The raw data entering this catalog 
has been reduced using the public software package grizli [33], which masks 
imaging artifacts, provides astrometric calibrations based on the Gaia DR3 
catalog, and shifts the images to a common pixel scale of 0."04/pixel using 
astrodrizzle. This catalog also includes the most recent updated photo- 
metric zero-points, and are all corrected for the Milky Way extinction. The 
catalog is published at this MAST DOI: 10.17909/g2av-£122. We use the 
aperture-matched photometry derived using a diameter of 0."5 from the cat- 
alog. A subset of the galaxies in our sample with NIRSpec prism spectra 
from the CEERS survey (CEERS-27789, -28612, -z8172, and -27453) are not 
covered by the JWST NIRCam mosaic. For those targets we use archival 
data from the Hubble Space Telescope (HST) processed in the same manner 
as the JWST mosaics described above, and published at this MAST DOI: 
10.17909/qhte-a822. 

The 1D spectra are extracted from a custom-made pipeline [34], that uti- 
lizes the Stage 2 output from the MAST archive. This code performs the 
standard wavelength, flat-field and photometric calibrations to the individual 
NIRSpec exposure files. We use the JWST reference files associated with the 
CRDS context jwst 1027.pmap. We generate the full combined 2D spectra 
and optimal [35] 1D extractions with scripts that extend the standard JWST 
pipeline functionality [34]. We further scale the overall flux densities of the 
spectra to the derived photometry, matching the integrated flux within the 
available passbands using a wavelength-dependent polynomial function. This 


ti.e., calibration levels 1 and 2 with jwst version 1.8.2 
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is to improve the absolute flux calibration of the spectra and take into account 
potential slit-losses, and is further used in the SED modelling of each source. 
An example is shown in Fig. 5 for the galaxy CEERS-z7832. The required 
scaling factors are typically of the order 1.2 — 2 and varies by 20-30% across 
the spectrum. The NIRSpec prism spectra typically cover the full wavelength 
range from 0.7m to 5.2um, except for RXJ-z8149 and CEERS-z7179, which 
have gaps at < 3um and 1.5—4ym, respectively. The spectral resolution range 
from R = 50 in the blue end to R = 400 in the red end [36]. The full set of 
reduced and photometrically-calibrated 1D NIRSpec/prism spectra are shown 
in Fig. 6. 


Lens models. To infer the relevant magnification factors for each source in 
the RXJ-2129 and Abell-2744 fields, we adopt the lens model presented by [37] 
for the RXJ-2129 galaxy cluster and consider the lens models compiled for the 
Hubble Frontier fields ? for the Abell 2744 cluster. We adopt the GLAFIC-v4 
lens model results in this work, but overall find that magnification factors are 
consistent within 5 — 1096 for the various online lens models. The magnification 
factors are summarized in Table 1, typically of the order u = 1 — 2 except for 
the galaxy RXJ-z9500 which has u = 19.2 + 3.6 [38]. 


SED modelling. The SEDs of each galaxy have all been modelled using 
BAGPIPES [10] to infer their physical properties such as the stellar masses and 
mass-weighted ages. T'his code uses the stellar population models from Bruzual 
& Charlot [39]. We fix the redshifts to zspec and include nebular emission via 
CLoupv [40]. We jointly model the photometric and spectroscopic data for 
each source, which greatly improves the statistical uncertainties from the pos- 
terior distributions. We further modify the default capabilities of BAGPIPES 
to incorporate a wavelength-dependent spectral resolution matching the NIR- 
Spec prism spectra. We allow the ionization parameter, U, to vary between 
—3 « log(U) « —1, to account for the typically higher ionization parame- 
ters at z = 6 [41]. We assume the attenuation curve from Salim et al. [42], 
which better represent low-mass galaxies, allowing the slope (6) to vary from 
—0.9 < ô < 0.1 and the bump strength (B) from 0 < B < 3. We assume Gaus- 
sian priors on the various dust parameters and infer Ay < 0.1 mag for the 
majority of galaxies. We assume a Kroupa 2001 initial mass function (IMF) 
[15] and set the metallicity to vary between 0 < Z/Ze < 0.5. 

Based on the accurate spectro-photometric modelling, we incorporate a 
more flexible, non-parametric SFH for each galaxy [43]. We find that this 


https: //archive.stsci.edu/prepds/frontier /lensmodels/ 
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effectively increases the inferred stellar masses and mass-weighted ages of the 
stellar population by 0.5-1 dex compared to more simple, constant SFHs, as 
also demonstrated in previous work [11, 44-46]. This model fits for quartiles of 
total mass formed, with a prior included to account for possible hidden mass 
due to “outshining” from the younger, burstier star-forming populations and 
includes an additional prolonged, less active period of star formation. Physi- 
cally, this corresponds to including the older, more massive stellar populations 
recovered from spatially-resolved analyses [12]. Table 1 reports the median and 
16th and 84th percentiles from the resulting posterior distributions on the stel- 
lar mass, dust attenuation, and mass-weighted ages of the stellar populations 
for each galaxy. 

We tried various sets of combinations using different parametrizations of 
the SFH and dust laws, including a constant and delayed SFH and a Salim 
[42] vs. Calzetti [47] dust attenuation curve. An exponentially declining SFH 
is likely not an accurate representation of these very high-redshift galaxies and 
may overestimate the inferred stellar masses with increasing SFRs. Moreover, 
the constant SFH is likely a too simple representation of the more complex 
SFHs expected for these galaxies, as predicted from simulations [48]. Assuming 
a Salim over a Calzetti-like attenuation curve effectively decreases the inferred 
Ay due to the steeper functional form of the Salim curve. The slope of the 
attenuation curve is related to the dust opacity and grain size distribution, 
which introduces and apparent dependence on the stellar mass as well [42] such 
that more more massive galaxies have shallower slopes. However, we found 
that our results were remarkably consistent overall (within 0.1 dex in stellar 


mass estimates), independent of the exact dust prescriptions used. 


Line flux measurements. To estimate the line fluxes for each galaxy, 
we model the most prominent redshifted nebular emission lines: H8, Hy, 
[O 111] 44959, 5007, and [O 11] 43727 imposed on the photometrically-calibrated 
spectra with simple Gaussian line profiles. The spectral resolution is insuffi- 
cient to determine more accurate line-profile kinematics. The redshift zspec and 
line widths are tied for all the lines, such that only the line fluxes can vary 
across transitions. We further fix the line flux ratio of the [O 111] A44960, 5008 
doublet to the theoretical predicted value of fiorm.-490/ ftormj-5oos = 1/3 in 
the fits, based on the relative transition probabilities from the excited "D, 
state to ?P, and ? P5, respectively [16]. 


Metallicity calibrations. To determine the gas-phase metallicity for each 
galaxy we use the strong-line diagnostics from Nakajima et al. [23], since we 
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are not able to resolve the auroral [O 111] 44363 line transition in all but one 
case. Instead, we chose the R3 diagnostic since it has a small empirical scatter 
(0.16 dex), and is supported by recent theoretical framework [23] and empiri- 
cal constraints based on direct metallicities estimates through the T.-method 
[24]. Further, the [O 1] A5008 and H8 emission lines are closely separated in 
wavelength space and this calibration therefore introduce less uncertainties 
due to the flux calibration and dust attenuation of the spectrum. Specifically, 
we adopt the “Large EW” calibration, representing the stronger ionization 
states of high-redshift galaxies [24], and also consistent with the inferred large 
H8 EWs > 200 À for the galaxies in our sample. The inferred metallicities are 
summarized in Table 1. Comparing the estimates from the R3 calibrations to 
the R2, R23, and O32 calibrations from [23], all requiring the additional detec- 
tion of [O r1] 43727, we find that they are overall consistent within 0.1-0.2 dex. 
For three of the galaxies in our sample, Abell2744-z7885, CEERS-z8612, and 
CEERS-z8172, the metallicity could not be constrained using the R3 diagnos- 
tic since their line ratios were at the turnover of this particular calibration. 
Instead we use a combination of the R2 and O32 calibrations since the derived 
line fluxes allow for more accurate metallicity estimates due to the distinct 
functional form of these diagnostics. 

For the one galaxy RXJ-29500, first presented by Williams et al. [38], 
we constrain the metallicity based on the marginal detection of the auroral 
[O 111] 44363 line at onn-me Z 5.0x 1071? ergs ^! cm™?, see Fig. 9. We follow 
the iterations outlined by Izotov et al. [49] and solve for the oxygen abundance 
as a function of electron temperature and density. We derive electron temper- 
atures T, X 2.8 x 10* K and gas-phase metallicities 12+ log(O/H) > 7.2 for all 
ne < 10* cm^?. This is consistent with the inferred estimate from the R3 cali- 
bration of 12 + log(O/H) = 7.56*015 and with previous results from Williams 
et al. [38]. 


Simulations of the mass-metallicity relation at z = 7— 10. In the 
Main Text, we found a large scatter in the predictions for the mass-metallicity 
relation of galaxies at z = 7 — 10 across a suite of state-of-the-art cosmological 
simulations. The predictions from the TNG50 simulations [extracted using 
the abundance ratios of gas cells weighted by the star formation rate within 
twice the half mass radius of galaxies; 22, 50] at z = 7 — 10 almost directly 
follow the observed mass-metallicity relation for galaxies at z — 0 [2], which 
is inconsistent with the overall evolution of this relation, even at z « 3. The 
FLARES simulations [measured by using the abundance ratios of gas particles 
weighted by their star formation rate within a 30 pkpc radius from the centre 
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of potential; 27, 51], when extrapolated to low stellar masses, seem to be 
marginally consistent with the observed mass-metallicity relation for galaxies 
at z = 7 — 10. However, the slope appears to be significantly steeper than 
105 Mo. 
It should be noted, however, that both the absolute normalization and the slope 


observed, effectively predicting more metal-enriched galaxies at M, = 
of the mass-metallicity relation are subject to the uncertainties in the metal 
yields of stars of all masses as well as the stellar feedback models implemented 
in the different simulations. For this reason, predicting the relative evolution of 
the mass-metallicity relation with redshift has been the most important so far. 
We also note that previous works, when comparing to theoretical predictions 
le.g. 28], assume that 35% of the total metal fraction is in oxygen when deriving 
the oxygen abundance (O/H), which will likely be different at these redshifts 
due to a—enhancement [52]. 

The Astraeus [25] and FirstLight [26] simulations appear to accurately 
match our empirical results. In the Astraeus simulation, masses of the metal 
and gas reservoirs are evolved by applying the semi-analytic galaxy evolution 
model described in [19] and [25] to the merger trees of halos identified in 
the dark matter-only Very Small MultiDark Planck N-body simulation. Gas- 
phase metallicities are then derived by taking the ratio between the oxygen 
and gas mass reservoirs, assuming the latest stellar metal yields from [53]. The 
FirstLight simulations output are taken directly from the literature [26], but 
corrected by 0.21 dex to match the most recent estimate of the solar metallicity, 
12 + log (O/H)c = 8.69 [54] (from the assumed 12 + log (O/H) = 8.9). 


Robustness of the high-redshift FMR offset. To compare the local 
FMR parametrized by Curti et al. [2] to the high-redshift observations, we 
converted their results based on the Chabrier [55] IMF to match the Kroupa 
IMF [15] assumed in this work. Effectively, this increases the stellar masses 
by z 10% and decreases the SFRs by ~ 20%. Since the functional form of 
the metallicity as a function of the stellar mass and SFR is prescribed as 
log(M,/Mo) — 0.55 x log(SFR/Moyr~') in the FMR, this change in IMF has 
a negligible effect. Even if the uncertainties associated with the metallicity 
predictions provided by the FMR parametrization in [2] increases to ~ 0.3 dex 
in the low-mass-high-SFR regime we are probing with this high-redshift galaxy 
sample, we still observe a systematic offset at 2c confidence. However, based 
on the high-redshift FMR predicted by the Astraeus simulations [25], the offset 
might increase even more at the low-mass-high-SFR regime, substantiating 


our results even more significantly. 
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As detailed above, assuming a simple constant SFH in the SED modelling 
of each galaxy will yield lower stellar mass estimates by 0.5-1 dex. Based on 
these estimates, which are effectively lower bounds, we still observe a factor 
of ~ 2 offset in the high-redshift galaxy sample from the local FMR at 20 
confidence. Accounting for a potential lower Ha-to-SFR calibration due to the 
higher efficiency of ionizing photon production in more metal-poor environ- 
ments [56], will systematically decrease the quantity log(M,/Mc) — 0.55 x 
log(SFR/Mayr—') by 0.2 dex. Even with these more conservative estimates 
we still recover a significant offset from local FMR in the metallicities of the 
high-redshift galaxy sample. 

Sanders et al. [5] established that the local FMR [1, 2] holds at z ~ 2.3 and 
z £4 3.3, based on an extensive sample of galaxies from the MOSDEF survey. 
'This suggests that galaxies evolve through a smooth baryonic growth process 
in a state of near-equilibrium from redshifts z — 0 — 3.3, i.e. from 2 Gyr after 
the Big Bang to the present. Comparing our observations to the recent sample 
of strong [O 11|-emitting galaxies at z = 5— 7 from Matthee et al. [8] in Fig. 10 
reveal no strong evolution from z % 6 to z ~ 8. The equilibrium state of galaxy 
growth and evolution thus appears to already be broken at z ~ 6, i.e. at 1 Gyr 
after the Big Bang. In Fig. 10, we also compare our observations to the three 
galaxies at z > 7 identified in the Early Science Observations (ERO), with 
inferred physical properties from [28]. These three sources all have metallicities 
derived from the direct 7T,-method, and are overall well in agreement with our 
results. The independent measurements at z — 5 — 7, and through the direct 
T.-method for the few galaxies at z > 7, thus greatly substantiate and support 
our results. Further, they imply that galaxies down to at least z ~ 6 are still 
subject to significant gas accretion from the intergalactic medium, causing a 
dilution of the chemical abundances in the interstellar medium of these sources. 
This is supported by the abundant neutral atomic gas reservoirs observed in 
galaxies at this epoch [30, 57]. The pristine, intergalactic gas is likely exhausted 
by z £ 3, about 1 Gyr later, at which point the local FMR is established. 
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Fig. 5: The original (blue) and rescaled (orange) 1D spectrum of the galaxy 
CEERS-27832. The NIRSpec prism spectrum has been scaled to match the 
NIRCam photometry to improve the accuracy of the flux calibration and 
account for potential slit-losses. The NIRCam filter curves used to photomet- 
rically calibrate the spectrum is shown by the grey passbands. The derived 
polynomial function, which we apply to correct the original spectrum, is shown 
in the top panel. 
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Fig. 6: Photometrically-calibrated 1D NIRSpec/prism spectra of the full 


sample of galaxies at z = 7 — 10. 
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Fig. 8: The photometrically-calibrated 1D spectrum of the galaxy CEERS- 
27832 (orange), imposed with the best-fit continuum and line emission from 
the BAGPIPES SED model (black). The grey region represents the uncertainty 
on the flux. The output posterior distributions from the best-fit SED on the 
stellar mass M,, mass-weighted age Tmass, dust attenuation Ay, and ionization 
parameter U, are shown at the top. 
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'Table 2: Line flux measurements. 


Galaxy ID Zspec HB [Om] A5007 [O 1] A3727 


RXJ-z9500 9.5008 18.01 + 1.61 91.44 + 1.90 4.94 + 1.87 
RXJ-z8149 8.1496 5.76 + 2.27 19.62 + 2.57 < 2.54 
RXJ-z8152 8.1523 9.37 + 1.90 53.99 + 2.23 < 5.21 
Abell-z7878 7.8783 < 8.01 28.37 + 4.58 9.67 + 4.55 
Abell-z7885 7.8854 7.13 X 2.27 77.69 + 2.67 8.74 + 2.66 
CEERS-z7105 7.1051 5.27 + 1.23 33.72 + 1.44 < 2.44 
CEERS-z7179 7.1792 7.84 + 0.73 36.38 + 0.84 


2.86 7.54 + 2.66 
1.56 2.73 + 1.39 
7.10 18.68 + 5.83 
1.23 4.33 + 1.05 
E 5.57 7.89 + 4.58 
1.86 3.81 + 1.62 
0.98 1.91 + 0.79 
1.35 3.50 + 1.12 
0.76 1.10 + 0.66 


CEERS-z7175 7.1751 14.4610 2.64 60.91 
CEERS-z7167 7.1683 3.01 + 1.38 19.73 
CEERS-z8684 8.6849 57.60 +5.95 461.03 
CEERS-z77T89 7.7805 « 2.61 19.99 
CEERS-z8718 8.7182 12.27 +4.68 62.04 
CEERS-z7832 7.8328 14.30 +1.60 87.42 
CEERS-z8612 8.6142 . 2.69 X 0.80 15.87 
CEERS-z8172 8.1724 8.94 + 1.12 64.09 
CEERS-z7453 7.4536 2.82 + 0.66 11.85 


Notes. The derived line fluxes are in units of 101° ergs™1 cm7? and have not been 
corrected for the magnification factor u. Upper limits denote 2e upper bounds on 


the line fluxes. 
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Fig. 9: (Top): Zoom-in on the photometrically-calibrated 1D spectrum of the 
galaxy RXJ-29500 (black), marking the most prominent emission line fea- 
tures. The typical strong nebular emission lines are detected, in addition to a 
marginal detection of the auroral [O 111] \4363 transition. The best-fit model 
of the continuum and line profiles are shown in red. (Bottom): The derived 
electron temperature T, and oxygen abundance 12 + log(O/H) for RXJ-z9500 
as a function of electron density ne, based on the direct T;-method [49]. 
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Fig. 10: The fundamental-metallicity relation (FMR) of galaxies, similar to 
Fig. 4. For comparison, we show the sample of [O1m]|-emitters at z = 5 — 7 in 
different mass bins [8], and the three galaxies at z > 7 with direct metallicity 


estimates from the T.-method [28]. 
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